% $Id: nlss.tex 204 2009-04-20 01:59:36Z jsibert $
%
% Author: David Fournier
% Copyright (c) 2008 Regents of the University of California
%

%\magnification=1400
%\def\mysection#1{{\noindent\bf\medskip #1\medskip}}
\def\wha{\widehat\alpha}
\mysection{Introduction}
An enduring characteristic of computer science and technology
appears to be that hardware develops more quickly than software,
so that it is good practice to devlopo software for tomorrows
hardware. Currently it appears that the next dvelopment in computer hardware will be the widesparead availability of mult-core processors so that very soon
processors with eight, sixteen, or more processing cores will be commonplace.
The obvvious software devlepment to take advantage of this hardware devlopment 
is parallel processing, that is computer programs which do there computaitons in parallel. However the devlopment of parallel programming software 
whihc could permit the user to write custom parllel programs the has been 
disappointing. There has been more success in developing parallel code for
particular ``generic'' applications such as linear algebra 
where the class of situations to be dealt  with is well understood in advance.
The lesson is that to extend the applcability of parallel processing 
one should pick a class of models or programming problems which is 
well understood and impelemtn the parallel aspects of the programming in such as
way that they they are as transparent as possible to the end user.
In this paper we describe such a programming paradigm for 
nonlinear (and linear) state spcae models. The programs themselves are
impelemtnbed using the random effects software module from
the AD Model Builder software package.

\def\Byt{\boldsymbol{y_t}}
mysection{The linear state space model}
To intorudce the ideas we begin with the linear state space model where the
calculations are very simple. We have adopted the  notation used in Harvey.
Let $\Byt$ be a multivariate time series

mysection{The nonlinear state space model}
The nonlinear state space model is specified by the set of equations
\begin{equation}
 \myeq{\alpha_t=g_t(\alpha_{t-1},\epsilon_t;\theta)}\label{nlss:xx1}
\end{equation}
where $1\le t\le n$ and $\epsilon_t$
is an n~dimensional  random vector with probability density function
$\Phi_t(\epsilon_t;\theta)$. 
The vectors $\alpha_t$ are the "states" the of system
which are not observed directly but affect the
probability distribution of a set of observations
$Y_t$ where  $Y_t$ is a random vector with
proability density function $F_t(y_t,\alpha_t;\theta)$ that is
\begin{equation}
 \myeq{Y_t \sim F_t(y_t,\alpha_t;\theta)}\label{nlss:xx1a}
\end{equation}
The $\theta$ are a relatively small numbder of parameters
which we want to estimate. For the moment consider them to be fixed.
Then $F_t(y_t,\alpha_t;\theta)$ is the probability density 
function for $Y_t$ if we know $\alpha_t$. However in the state space model
the $\alpha_t$ are unknown. All we know is that the $\alpha_t$ evolve 
according to \ref{nlss:xx1}.
 To calculate the probability density
function $L(y_t,;\theta)$ for the $Y_t$ given only the 
 relationship \ref{nlss:xx1}
it is necessary to integrate over all the $\epsilon$
\begin{equation}
 \myeq{L(y_t;\theta)=\int \prod_t F_t(y_t,\alpha_t(\epsilon);\theta)\prod_t \Phi_t(\epsilon_t;\theta) d\epsilon}
\end{equation}
Note that in general $\alpha_t$ depends on all the $\epsilon_{t^\prime}$
for $t^\prime\le t$ so that this integral is difficult to evaluate. 
We shall construct a new parameterization with respect to which the
integral will have a simpler (but still difficult) form. 

%However first 
%we will derive the equations for the Kalman filter for linear state space models and show how that connects via optimziation over $\alpha$ to our
%general approach.
%\mysection{The Kalman filter for linear state space models}
%\begin{equation}
% \myeq{
%  (y-Z\alpha-d)^\prime H^{-1}(y-Z\alpha-d)+(\alpha-e)^\prime S^{-1}(\alpha-e)
%}\label{nlss:xx100}
%\end{equation}
%Consider this as a quadratic function in $\alpha$.
%\begin{equation}
% \myeq{
%  h(\alpha)=A+ B\alpha +{1\over2}\alpha^{\prime}C\alpha 
%}\label{nlss:xx101}
%\end{equation}
%Now let $\wha=-C^{-1}B$, so that
%\begin{equation}
% \myeq{
%  h(\wha+\delta)=A+ B\wha +{1\over2}\wha^{\prime}C\wha 
%     +{1\over2}\delta^{\prime}C\delta
%}\label{nlss:xx102}
%\end{equation}
%which simplifies to
%\begin{equation}
% \myeq{
%  h(\wha+\delta)=A -{1\over2}\wha^{\prime}C\wha 
%     +{1\over2}\delta^{\prime}C\delta
%}\label{nlss:xx102}
%\end{equation}
%so the integral becomes
%\begin{equation}
% \myeq{
%   \int \exp( -A +{1\over2}\wha^{\prime}C\wha) 
%     \exp(-{1\over2}\delta^{\prime}C\delta)\,d\delta
%}\label{nlss:xx103}
%\end{equation}
%or
%\begin{equation}
% \myeq{
%   |C|^{-{1\over2}}\exp( -A +{1\over2}\wha^{\prime}C\wha) 
%}\label{nlss:xx103}
%\end{equation}
%since $\wha=-C^{-1}B$ this simplifies to 
%\begin{equation}
% \myeq{
%   |C|^{-{1\over2}}\exp( -A +{1\over2}B^{\prime}C^{-1}B) 
%}\label{nlss:xx103}
%\end{equation}
%solving for $A$ and $B$ we get
%\begin{equation}
% \myeq{
%   A=(y-d)^\prime H^{-1} (y-d) +e^\prime S^{-1} e
%}\label{nlss:xx103}
%\end{equation}
%and
%\begin{equation}
% \myeq{
% B=(y-d)^\prime H^{-1}Z+e^\prime S^{-1}
%}\label{nlss:xx104}
%\end{equation}
%and
%\begin{equation}
% \myeq{
% C=Z^\prime H^{-1}Z+S^{-1}
%}\label{nlss:xx105}
%\end{equation}
%and
%\begin{equation}
% \myeq{
% B^\prime C^{-1}B=\Big((y-d)^\prime H^{-1}Z+e^\prime S^{-1}\Big)\Big\{Z^\prime H^{-1}Z+S^{-1}\Big\}^{-1}
%  \Big(Z^\prime H^{-1}(y-d) +S^{-1}e\Big)
%}\label{nlss:xx106}
%\end{equation}
%


Now for any value of the parameters $\alpha_{t-1}$
assume that it is possible to solve the  equations\ref{nlss:xx1}
for  $\epsilon_t$ in terms of $\alpha_t$ ie
\begin{equation}
 \myeq{\epsilon_t=h(\alpha_t,\alpha_{t-1})}\label{nlss:xx3}
\end{equation}
so that for any fixed value of $\alpha_{t-1}$ the map \ref{nlss:xx3}
is a diffeomorphism. The Jacobian matrix $J(\alpha)$ of this map
is a banded upper triangular matrix of the form
\[
\left | \begin{array}{ccccc}
   D_1h_1(\alpha_1,\alpha_0) & D_2h_2(\alpha_2,\alpha_1)  & 0 &  & 0 \\
     0                    & D_1h_2(\alpha_2,\alpha_1)  &  D_2h_3(\alpha_3,\alpha_2) & & 0  \\ 
     0                    &  0                        &  D_1h_3(\alpha_3,\alpha_2)  & & 0\\
                          &                           &                & \ddots \\
     0                    & 0                         &    0           &           & D_2h_n(\alpha_n,\alpha_{n-1}) \\
     0                    & 0                         &    0           &           & D_1h_n(\alpha_n,\alpha_{n-1}) 
 
   \end{array}
  \right |
\]
where $D_1h_t(\alpha_t,\alpha_{t-1})$ is the (square) matrix of partial derivatvies of $h_t$ with respect to  $\alpha_t$ and 
 $D_2h_t(\alpha_t,\alpha_{t-1})$ is the matrix of partial derivatives 
of $h_t$ with respect to  $\alpha_{t-1}$. 
Let $\Phi_t(\epsilon_t)$ be the proabability density function of
the random vector $\epsilon_t$.  We can use the $h_t$  to
repararameterize the integral  in terms of the $\alpha_t$. 
This reparameterization involves the determinant of the Jacobian matrix, $J$.
 
\begin{equation}
 \myeq{
    \int \prod_t F_t(y_t,\alpha_t)
    |J(\alpha)|^{-1}\prod_t \Phi_t(h_t\big(\alpha_t,\alpha_{t-1})\big) d\alpha}\label{nlss:xx5}
\end{equation}
\mysection{Approximate integration via the Laplace approximation}
  The laplace approixmation for the integral \ref{nlss:xx5}
  involves approximating the integrand by the second order Taylor expansion 
  in $\alpha$ evaluated at the maximzing value $\widehat \alpha$.
    It is more convenient to maximize the log of the integrand or
   equivalently to minimize minus the log of the integrand.
    Let $f_t(y_t,\alpha_t)=\ln\big(F_t(y_t,\alpha_t)\big)$
    and $\phi_t(h_t\big(\alpha_t,\alpha_{t-1})=
            \ln\big(\Phi_t(h_t\big(\alpha_t,\alpha_{t-1})\big)$
We can write the integral as
\begin{equation}
 \myeq{
    \int \exp\Big( -\left\{-\sum_t f_t(y_t,\alpha_t)
    +\ln\big(|J(\alpha)|\big)
     -\sum_t \phi_t(h_t\big(\alpha_t,\alpha_{t-1})\big) \right\}\Big)
    \,d\alpha}
    \label{nlss:xx6}
\end{equation}
The computationally difficult term in \ref{nlss:xx6} is 
    $\ln\big(|J(\alpha)|\big)$.
Since $J(\alpha)$ is upper diagonal this is given by 
\begin{equation}
 \myeq{
    \ln\big(|J(\alpha)|\big)=\sum_t \ln(|D_1h_t(\alpha_t,\alpha_{t-1}|)}
    \label{nlss:xx7}
\end{equation}
which only involes the diagonal terms so that is it is not necessary to 
calculate the terms $|D_2h_t(\alpha_t,\alpha_{t-1})|$.
Gathering together terms with the sam $t$ subscript we obtain
\begin{equation}
 \myeq{
    \int \exp\Big( -\left\{-\sum_t f_t(y_t,\alpha_t)
     +\ln(|D_1h_t(\alpha_t,\alpha_{t-1})|
     -\phi_t(h_t\big(\alpha_t,\alpha_{t-1})\big) \right\}\big)d\alpha}
    \label{nlss:xx8}
\end{equation}
Let $\psi_t$ be defined by
\begin{equation}
 \myeq{
   \psi_t(\alpha_t,\alpha_{t-1};\theta)=-f_t(y_t,\alpha_t;\theta)
     +\ln\big (|D_1h_t(\alpha_t,\alpha_{t-1};\theta)|\big)
     -\phi_t\big (h_t(\alpha_t,\alpha_{t-1};\theta)\big)}
    \label{nlss:xx9}
\end{equation}
and rewrite the integral in terms of the $\psi_t$
\begin{equation}
 \myeq{
    \int 
      \exp\big( -\bigg\{-\sum_t \psi_t(\alpha_t,\alpha_{t-1};\theta)\bigg\}\big)
    \,d\alpha}
    \label{nlss:xx10}
\end{equation}
Now setting 
\begin{equation}
 \myeq{
   \Psi(\alpha;\theta)=\sum_t \psi_t(\alpha_t,\alpha_{t-1};\theta)}
    \label{nlss:xx11}
\end{equation}
so that the integral becomes
\begin{equation}
 \myeq{
    \int \exp\big( -\Big\{-\Psi(\alpha;\theta)\Big\}\big)\,d\alpha}
    \label{nlss:xx12}
\end{equation}
Let $\widehat\alpha$ be the value of $\alpha$ which minimizes
    $\Psi(\alpha;\theta)$ and approximate $\Psi(\alpha;\theta)$ 
 by its second order Taylor expansion at 
$\widehat\alpha$ we obtain
\begin{equation}
 \myeq{
    \int \exp\big( -\Big\{-\Psi(\widehat\alpha;\theta) 
     - D^2_\alpha\Psi(\wha;\theta)\Big\}\big)\,d\alpha}
    \label{nlss:xx13}
\end{equation}
which is equal to $\exp(\Psi(\wha;\theta)|-D^2_\alpha\Psi(\wha;\theta)|^{1/2}$
where $|-D^2\Psi(\wha;\theta)|$ is the 
determinant of the Hessian of $-\Psi(\wha;\theta)$.

$\Psi$ is separable function of $\alpha$that is it is
the sum of $n$ terms each of which involves at most two of the $\alpha_t$ and
each $\alpha_t$. This implies that the Hessian 
matrix is  banded (actually block diagonal). The banded structure can be 
exploited in two ways. Calculating the determinant of a block diagonal matrix
is very efficient. Also solving the system of equations 
$u=H^{-1}v$ where $H$ is a positive definite symmetric block diagonal matrix
can be done very efficiently. This latter property can be exploited to
help solve for $\wha$ via the newton-raphon method and its variants.
 
These techniques are incorporated into the AD Model Builder
random effects software module. It is only necessary for the user to
formulate the problem in the correct fashion using the 
{\tt SEPARABLE\_FUNCTION} option. 

\beginexampledf
DATA_SECTION
  init_int n
  init_matrix Y(0,n,1,3)
PARAMETER_SECTION
  init_vector alpha0(1,3)
  init_bounded_number log_sigma_e(-2.0,5.0)
  init_bounded_number log_sigma_y(-2.0,5.0)
  random_effects_matrix alpha(1,n,1,3)
  matrix eps(1,n,1,3)
  objective_function_value f
PROCEDURE_SECTION
  f0(alpha0,log_sigma_y);
  int i;
  //fx(alpha(1),alpha0);
  f1(alpha(1),alpha0,log_sigma_y,log_sigma_e,1);
  for (i=2;i<=n;i++)
  {
    f1(alpha(i),alpha(i-1),log_sigma_y,log_sigma_e,i);
  }

SEPARABLE_FUNCTION void f0(const dvar_vector& alpha0,const prevariable& log_sigma_y)
  dvariable s_y=exp(log_sigma_y);
  f+=3.0*log_sigma_y+0.5*norm2((Y(0)-alpha0)/s_y);
  
SEPARABLE_FUNCTION void f1(const dvar_vector& alphai,const dvar_vector& alphai1,const prevariable& log_sigma_y,const prevariable log_sigma_e,int i)
  dvariable s_e=exp(log_sigma_e);
  dvariable s_y=exp(log_sigma_y);
  dvar_vector epsi=alphai-alphai1;
  f+=3.0*log_sigma_e+0.5*norm2(epsi/s_e);
  f+=3.0*log_sigma_y+0.5*norm2((Y(i)-alphai)/s_y);

\endexampledf

\mysection{Advantages of the holistic approach}

For the linear model with normally distributed random variables
these are the full information maximum likelihood estimates. 

For nonlinear models and/or with non-normal random variables
 the holistic approach produces better estimators than the Kalman filter.

The holistic approach is easily parallelizable so it 
is a  good candidate for a computational paradigm which can exploit
future multi-core computers.

Consider expression \ref{nlss:xx10}. 
To calculate the Laplace approximaiton it 
is necessary to maximize this expression with respect to the parameters 
$\alpha_t$. This is the inner optimization. It has a separable
strcture that is it is a sum of  summands $\psi_t$
where each summand depends only on
$\alpha_t$, $\alpha_{t-1}$, and
$\theta$. As functions of these parameters the calcuations for each $\psi_t$
are independent of each other and so can be trivially parallelized,
and depend on only a relatively small number of parameters.
After (and optionally during) the inner optimziation it is 
also necessary to calculate the banded Hessian matrix, and to
caluclate its determinant. These calculations are also 
parallelizable and constitute all the calculations necessary tocarry out 
one calculation  of $L(\theta)$. Finally we need to calculate the
derivatives of $L(\theta)$ with respect to the parameters $\theta$.
These derivatives can also be calculated within the 
parallelizing structure used for the inner optimziation,
and share the reduction in computational complexity due to the separable
structure in the model.

\def\Nipjp{N_{i+1,j+1}}
\def\Nipm{N_{i+1,m}}
\def\Nij{N_{ij}}
\def\nij{n_{ij}}
\def\Ei{E_i}
\def\Nio{N_{i1}}
\def\Nim{N_{im}}
\def\Nnpj{N_{n+1,j}}
\def\nnpj{n_{n+1,j}}
\def\Nimm{N_{i,m-1}}
\def\Zij{Z_{ij}}
\def\Fij{F_{ij}}
\def\Zim{Z_{im}}
\def\Zimm{Z_{i,m-1}}
\def\qip{q_{i+1}}
\def\qi{q_i}
\def\ui{u_i}
\def\sj{s_j}
\def\deltai{\delta_i}
\def\di{d_i}
\def\epsij{\epsilon_{ij}}
\def\gami{\gam_i}
\def\lami{\lambda_i}
\def\etai{\eta_i}
\def\Cijobs{C^{obs}_{ij}}
\def\Cij{C_{ij}}

\mysection{A simple age-structured model}
This is an example of a simple model used in fisheries
stock assessment. There are $n$ years of catch data and the 
population is divided up into $m$ age classes where the last age 
class consists of all individual of age $m+a_0$ and older, 
where $a_0$ is the age of the
individuals in the first age class.
First we formulate the model in a manner which follows the 
stanbdard state space approach. 
Let $\Nij$ for $1\le i\le n+1$ and $1\le j\le m$ be the number of age class 
N indiviudual at the beginning of year $i$. 
For $1\le i\le n$ and $1\le j\le m$ the instantaneous total mortality rate
$\Zij$ is related to the $\Nij$ by the relationships
\begin{equation}
 \myeq{
   \Nipjp=\Nij\exp(-\Zij) \qquad \hbox{for}\qquad 1\le i\le n,\quad 1\le j<m-1}
    \label{nlss:xx200}
\end{equation}
and since the last age class is all the older fish
\begin{equation}
 \myeq{
   \Nipm= \Nimm\exp(-\Zimm)+ \Nim\exp(-\Zim)\qquad \hbox{for}\qquad 1\le i\le n}
    \label{nlss:xx201}
\end{equation}
Now the conventional way to parameterize the model is to start with
the number of fish at the beginning of year $1$ 
and the number of age class 1 fish at the beginning of each
$\Nio$ (recruitment) 
and to use the 
relationships \ref{nlss:xx200} and \ref{nlss:xx201} to compute
the other $\Nij$. 
It is of course necessary to sup[[ly values for the $Zij$, which we will deal 
with below, but the point is that one
starts with the fish population in a cetain "state" at the beginning of year
$1$ and calculates the later states as a function of this initial state and
its transitions from one year to the next. Now assume that
the total mortltiy can be decomposed into that due to fishing, $\Fij$
and that due to all other factors, $M$. 
where $\Fij$ is referred to as 
the instantaneous fishing mortality rate
and $M$ is referred to as the
the instantaneous natural mortality rate.
It is assumed that
\begin{equation}
 \myeq{
   \Zij=\Fij+M 
   }
    \label{nlss:xx202}
\end{equation}
It is further assumed that
\begin{equation}
 \myeq{
   \Fij = \qi \sj \di \Ei \exp(\epsij) 
   }
    \label{nlss:xx203}
\end{equation}
where the $\Ei$ are exogenous parameters referred to as the fishing effort.
\begin{equation}
 \myeq{
   \qip=\qi\exp(\deltai) \qquad \hbox{for} \quad 1<i<=n
  }
    \label{nlss:xx204}
\end{equation}
In addition assume the the recruitments $\Nio$ satify the relationship 
\begin{equation}
 \myeq{
   \Nio=R\exp(\lami) \qquad \hbox{for} \quad 1<i<=n
  }
\end{equation}

\begin{equation}
 \myeq{
   \di=\exp(\etai) \qquad \hbox{for} \quad 1<i<=n
  }
    \label{nlss:xx205}
\end{equation}
The parameters $\etai$, $\epsij$, $\deltai$, 
and $\lami$ are the random effects, 
that is, They are the random variables which we wish to integrate over to
obtain the marginal distribution.

Let $\Cij$ be the catch of age class $j$ fish in year $i$.
which are given by the relationship.
\begin{equation}
 \myeq{
   \Cij={\Fij\over \Zij}\big(1-\exp(-\Zij)\big)\Nij
     \qquad \hbox{for} \quad 1\le i\le n \quad 1\le j\le m
  }
  \label{nlss:xx206}
\end{equation}
We also need some observational data to fit the model. We shall assume that
we have estimates of the catch at age $\Cijobs$.

\mysection{Wholistic parameterization for the age-structured model}
The total number of random variables is
$nr$ $\epsij$ parameters, $n$ $\etai$ parameters,
$n$ $\lami$ parameters, and $n-1$ $\deltai$ parameters,
so there are $nr+3n-1$ random variables in total.

We shall reparameterize the random variables by 
$N_{11}$, 
$\Nij$ for $2\le i\le n+1$, $1\le j\le m$, 
$\Nnpj$ for  $2\le j\le m$, 
$d_i$ for $1\le i\le n$,
$\qi$ for $2\le i\le n$, and $\ui$ for $1\le i\le n$.
where the $\ui$ are needed to split up the $m$'th age class.



